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Introduction 


Motivation 


This course focuses on algorithms in signal processing and 
communications. When training data is available in such systems, we can 
process the data by first training on historical data and then running a 
Bayesian scheme, which relies on the statistics being known. A similar 
Bayesian approach can also be used when the statistics are approximately 
known. 


For example, consider lossy image compression. It is well known that 
wavelet coefficients of images have a distribution that resembles Laplace, 
Equation: 


i(2) = cyce lel 


and the coefficients are approximately independent and identically 
distributed (i.i.d.). A well-known approach to lossy image compression is to 
first compute the wavelet coefficients and then compress them using a lossy 
compressor that is designed for Laplace i.i.d. coefficients [link]. In this 
approach, the training consists of the body of art that realizes that wavelet 
coefficients are approximately Laplace i.i.d., and the Bayesan algorithm is a 
lossy compressor that is designed for this distribution. 


However, sometimes the statistics of the input (often called the source) are 
completely unknown, there is no training data, or there is great uncertainty 
in the statistics. For example, in lossless data compression, we do not know 
whether a file is an executable, source code, a DNA sequence, contains 
financial transactions, is text, etc. And even if we know that the file 
contains text, it has been noted that even different chapters that appear in 
the same book that is written by the same author may contain different 
Statistics. 


For this latter set of problems, the Bayesian approach is useless, because 
there is no training data. A good alternative approach to Bayesian 


algorithms is to use universal algorithms [link], [link]. These algorithms 
have good performance irrespective of statistics. In fact, in some cases, 
these algorithms can achieve (with equality) the theoretically optimum 
performance in an appropriate asymptotic framework. 


In lossless compression, universal algorithms have had great impact. For 
example, the Lempel-Ziv family of algorithms [link], [link] asymptotically 
achieve the entropy rate [link], [link], which is the best possible 
compression ratio achievable in lossless compression, despite not knowing 
the input statistics. Additionally, the Lempel-Ziv algorithms allow efficient 
implementation. 


Overview 


The goal of this course is to study universal algorithms, starting from the 
well-trodden material in lossless compression, and later discussing 
universal algorithms in other areas of communications and signal 
processing. Let us overview the material studied during the course. We 
begin in [link] with a review of some information theory material, including 
typical sequences and source coding, in order to provide sufficient 
background. Next, [link] statistical models for data will be described. [link] 
then presents techniques for universal lossless compression of parametric 
sources. One approach to universal compression of parametric sources is 
minimum description length, which compresses the data in order to 
minimize for the sum of the complexity of the model and the complexity of 
the data given the parameters of the model. Minimum description length 
has been used with context tree models to provide universal contextual 
prediction; context tree approaches are detailed in [link]. We then switch 
gears in [link] and move beyond lossless compression; universal lossy 
compression and signal reconstruction are described in detail. Finally, [link] 
describes Lempel-Ziv algorithms for universal lossless compression based 
on parsing an input sequence. For convenienve, notation is summarized in 
[link]. 


This manuscript is a work in progress, and we expect to expand and 
improve it during future teachings of the course. 


Background 


Convergence of random variables 


We include some background material for the course. Let us recall some 
notions of convergence of random variables (RV's). 


¢ A sequence of RV's {2n},,51 converges in probability if 
Ve > 0, lim,_,..supPr (|r, — %|> €) = 0. We denote this by 
P. 
Li 
¢ A sequence of RV's {2n},,5, converges to % with probability 1 if 
w.p.1 
Pr {1, X2,... ‘limp. 2n = £} = 1. We denote this by x, —> & 
ore, = aos 
e A sequence of RV's {Znhn>1 converges to % in the te sense if 


lp 
E||z, — %|?| — 0. We denote this by x, — &. 


M.S. 
For example, for p = 2 we have mean square convergence, 2, —> Z. For 
p= 2, 
Equation: 
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Elzn, —,z|P" = E(|zn, ie, < (Elzn _ z|”) 


Ly tat 
Therefore, x, — % yields x, —> 2. Note that for convergence in @, sense, 
we have 
Equation: 


E\z, — =| 
ee ee 


Pr (2, — &|> €) — 0. 
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Typical Sequences 


The following material appears in most textbooks on information theory 
(c.f., Cover and Thomas [link] and references therein). We include the 
highlights in order to make these notes self contained, but skip some details 
and the proofs. Consider a sequence x = x” = (21, 2%2,...,2n), where 

x; € @, a is the alphabet, and the cardinality of a is r, ie., |a| = r. 


Definition 1 The type of x consists of the empirical probabilities of 
symbols in z, 
Equation: 


where 7, (a) is the empirical symbol count, which is the the number of 
times that a € @ appears in z. 


Definition 2 The set of all possible types is defined as P,,. 


Example: 
For an alphabet a = {0, 1} we have 


P, = {(2, *), (4, =),..., (4, ©) }. In this case, 


non i? n? 


P,|=n-+1. 


Definition 3 A type class T, contains all x’ € a”, such that P,=P,,, 
Equation: 


Lea yVaie 6a i PvP): 


Example: 
Consider a = 1, 2,3 and x = 11321. We have n = 5 and the empirical 


counts are nz = (3,1, 1). Therefore, the type is P, = (2, = +) and the 


type class T’, contains all length-5 sequences with 3 ones, 1 two, and 1 
three. That is, J, = {11123, 11132, ..., 32111}. It is easy to see that 


IT |= yen = 20. 


Theorem 1 The cardinality of the set of all types satisfies 
\Pal< (m +1)". 


The proof is simple, and was given in class. We note in passing that this 
bound is loose, but it is good enough for our discussion. 


Next, consider an i.i.d. source with the following prior, 
Equation: 


Q(«) =] @ (a). 


i=1 


We note in passing that i.i.d. sources are sometimes called memoryless. Let 
the entropy be 
Equation: 


H (Px) = ~Zaca 2 tog (22 ) | 


where we use base-two logarithms throughout. We are studying the entropy 
H(P,,) in order to show that it is the fundamental performance limit in 
lossless compression. »/ find me 


We also define the divergence as 
Equation: 


DiPs | Q:) = XacaPs log (=). 


x 


It is well known that the divergence is non-negative, 
Equation: 


D(Pz || Qz) > 0. 


Moreover, D(P || Q) = 0 only if the distributions are identical. 


Claim 1 The following relation holds, 
Equation: 


Q (x) = 2-H (Ps) + D(Pel@@))] 


The derivation is straightforward, 
Equation: 


Q(x) = TreaQ(a)™ 
y) SacaNe (a)logQ(a) 
n5'P,(a) (log$-HogP) 


—n[H(Px)+D(P2||Q(x))] 


2 
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Seeing that the divergence is non-negative [link], and zero only if the 
distributions are equal, we have Q (x) < P, (x). When P, = Q the 
divergence between them is zero, and we have that 

P, (2) = Qz = 2-7H (Ps), 


The proof of the following theorem was discussed in class. 


Theorem 2 The cardinality of the type class T'(P,) obeys, 
Equation: 


(n a i ; gn (Pz) ra IT (Pz)| < gnH(P,) 


Having computed the probability of x and cardinality of its type class, we 
can easily compute the probability of the type class. 


Claim 2 The probability Q(T (P,) of the type class TP.) obeys, 
Equation: 


(n ab 1") : g-nD(P2||Qz) < Q (T (P,)) a gD Pr||Qz) | 


Consider now an event A that is a union over T(P,,). Suppose T(Q) Z A, 
then A is rare with respect to (w.r.t) the prior Q. and we have 

limn+oo Q (A) = 0. That is, the probability is concentrated around Q. In 
general, the probability assigned by the prior @ to an event A satisfies 
Equation: 


Q(A) = YeesQ (x) = “7(P,)cAQ (T (P:)) 
Dryp,yca2 MO Pai@ 
go r-minpeaD(PIQ), 


where we denote a,, = b, when log (=) —+ 0. 


Fixed and Variable Length Coding 


Fixed to fixed length source coding: As before, we have a sequence x of 
length n, and each element of x is from the alphabet a. A source code maps 
the input x” € r” toa set of 22” bit vectors, each of length Rn. The rate R 
quantifies the number of output bits of the code per input element of z. 
[footnote] That is, the output of the code consists of nF bits. If n and R is 
fixed, then we call this a fixed to fixed length source code. 

We assume without loss of generality that Rn € Z. If not, then we can 
round Rn up to | Rn], where |-| denotes rounding up. 


The decoder processes the nR bits and yields & € a”. Ideally we have that 
£ = x, but if 2”” < r” then there are inputs that are not mapped to any 


output, and & may differ from x. Therefore, we want Pr (@ i x) to be 
small. If R is too small, then the error probability will go to 1. On the other 
hand, sufficiently large R will drive this error probability to 0 as n is 
increased. 


If log (r) > Rand Pr (@ f x) is vanishing as 7n is increased, then we are 
compressing, because Qlog(r)n — pr s+ QR. where r” is the number of 
possible inputs a and there are 2%” possible outputs. 


What is a good fixed to fixed length source code? One option is to map 
2” _ 1 outputs to inputs with high probabilities, and the last output can be 
mapped to a “don't care" input. We will discuss the performance of this 
style of code. 


An input x € r” is called 5-typical if Q (x) > 2~(#+9)", We denote the set 
of 6-typical inputs by Tg (6), this set includes the type classes whose 
empirical probabilities are equal (or closest) to the true prior Q(x). Note 
that for each type class T’,, all inputs x’ € T, in the type class have the 
same probability, ie, Q (x’) = Q (x). Therefore, the set Tg (6) is a union 
of type classes, and can be thought of as an event A ((link]) that contains 
type classes consisting of high-probability sequences. It is easily seen that 
the event A contains the true i.i.d. distribution Q, because sequences whose 
empirical probabilities satisfy P, = Q also satisfy 

Equation: 


Using the principles discussed in [link], it is readily seen that the probability 
under the prior Q of the inputs in Tg (4) satisfies Q(T, (6)) = Q(A) > 1 
when n — oo. Therefore, a code @ that enumerates Tg (6) will encode z 
correctly with high probability. 


The key question is the size of @, or the cardinality of Tg (6). Because 
each z € Tg (6) satisfies Q (x) > 2(-#+4)", and d2eT (8) Q (x) < 1, we 
have Te (5)| < 2\7+9)”, Therefore, arate R > H + 6 allows near- 


lossless coding, because the probability of error vanishes (recall that 
Q G (5))°) —+ 0, where (-)© denotes the complement). 


On the other hand, arate R < H — 6 will not allow lossless coding, and the 
probability of error will go to 1. We will see this intuitively. Because the 
type class whose empirical probability is Q dominates, a type class T, 
whose sequences have larger probability, e.g., Q (x) > 2-4-9)" will have 
small probability in aggregate. That is, 
Equation: 
Nn—-0o 
Q (x) — 0. 
r:Q(x)>2-"(4—9) 


In words, choosing a code @ with rate R = H — 6 that contains the words 
x with highest probability will fail, it will not cover enough probabilistic 
mass. We conclude that near-lossless coding is possible at rates above H 
and impossible below H. 


To see things from a more intuitive angle, consider the definition of entropy, 
H (Q) = — daca @ (a) log (Q (a)). If we consider each bit as reducing 
uncertainty by a factor of 2, then the average log-likelihood of a length-n 
input x generated by Q satisfies 

Equation: 


E\— log (Pr (x))]| = z|- log (TI Pr @) 
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Because the expected log-likelihood of x is nF, it will take nf bits to 
reduce the uncertainty by this factor. 


Fixed to variable length source coding: The near-lossless coding above 
relies on enumerating a collection of high-probability codewords Tg (6). 
However, this approach suffers from a troubling failure for z ¢ Tg (6). To 
solve this problem, we incorporate a code that maps z to an output 
consisting of a variable number of bits. That is, the length of the code will 
be approximately nH on average, but could be greater or lesser. 


One possible variable length code is due to Shannon. Consider all possible 
x € a”. For each z, allocate |— log (Q())| bits to x. It can be shown that 
it is possible to construct an invertible (uniquely decodable) code as long as 
the length of the code /(z) in bits allocated to each z satisfies 


Equation: 
ye aE 
x 


This result is known as the Kraft Inequality. Seeing that 


Equation: 
Soe = Seven 
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we see that the length allocation we suggested satisfies the Kraft Inequality. 
Therefore, it is possible to construct an invertible (and hence lossless) code 
with lengths upper bounded by 

Equation: 


ly = |— log (Q(a))] < — log (Q(a)) +1, 


and we have 
Equation: 


E{l(z)| < E[— log (Q(z))| +1 =nH +1. 


This simple construction approaches the entropy up to 1 bit. 


Unfortunately, a Shannon code is impractical, because it requires to 
construct a code book of exponential size |a|". Instead, arithmetic 

codes [link] are used; we discussed arithmetic codes in detail in class, but 
they appear in all standard text books and so we do not describe them here. 


Source models 


For iid. sources, D(P; (x”) ||P2 (z") ) = nD(P; (2;) ||-P2 (x:)), which means that the 
divergence increases linearly with n. Not only does the divergence increase, but it does so by a 
constant per symbol. Therefore, based on typical sequence concepts that we have seen, for an 
x” generated by P, its probability under P2 vanishes. However, we can construct a distribution 
Q whose divergence with both P; abd P» is small, 


Equation: 
n 1 n 1 n 
ie = 5 Pile )+ 3 Pa (x Ve 

We now have for P;, 
Equation: 

1 ra 1 P, x” 

2 p(pp|i@) = 2#]Iog ——2@) 

n nm |” EP, (o) + }P2(2”) 
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On the other hand, + D(P; (x?)||Q (27)) > 0 [link], and so 
Equation: 


= > —D(P, (23) |IQ (22) > 0. 


By symmetry, we see that Q is also close to P, in the divergence sense. 


Intuitively, it might seem peculiar that @ is close to both P; and P, but they are far away from 
each other (in divergence terms). This intuition stems from the triangle inequality, which holds 
for all metrics. The contradiction is resolved by realizing that the divergence is not a metric, and 
it does not satisfy the triangle inequality. 


Note also that for two i.i.d. distributions P; and Ps, the divergence 
Equation: 


D(P, (#") || Ps (2")) = nD(Pr ||P2) 


is linear in n. If Q were i.i.d., then D(P; (x”) || Q (x})) must also be linear in n. But the 
divergence is not increasing linearly in n, it is upper bounded by 1. Therefore, we conclude that 
Q(-) is not an iid. distribution. Instead, Q is a distribution that contains memory, and there is 
dependence in Q between collections of different symbols of x in the sense that they are either 
all drawn from P, or all drawn from P 2. To take this one step further, consider K sources with 
Equation: 


then in an analogous manner to before it can be shown that 
Equation: 


D(P, (#)|I@ (#8) < — log (K). 


Sources with memory: Instead of the memoryless (i.i.d.) source, 
Equation: 


let us now put forward a statistical model with memory, 
Equation: 


Stationary source: To understand the notion of a stationary source, consider an infinite stream 
of symbols, ...,£_1, 20, £1, .... A complete probabilistic description of a stationary distribution 
is given by the collection of all marginal distribution of the following form for all ¢ and n, 
Equation: 


oe Oe, Cee (xt, Ltt1ys+ »Lt4+n-1)- 


For a stationary source, this distribution is independent of t. 


Entropy rate: We have defined the first order entropy of an i.i.d. random variable [link], and let 
us discuss more advanced concepts for sources with memory. Such definitions appear in many 
standard textbooks, for example that by Gallager [link]. 


1. The order-n entropy is defined, 
Equation: 


1 1 
A, = et (ety ++ +9 Bn) a ~— Ellog (P(1,---,2n))]- 


2. The entropy rate is the limit of order-n entropy, H =lim,_,.. H,,. The existence of this 
limit will be shown soon. 


3. Conditional entropy is defined similarly to entropy as the expectation of the log of the 
conditional probability, 


Equation: 
1 
A Nhe Ve iz Ps sal oa Ge ceny  eeemremare soa 
where expectation is taken over the joint probability space, P(x1,..., Zn). 
The entropy rate also satisfies H =lim,_,.. H (ia isa). 


Theorem 3 For a stationary source with bounded first order entropy, H, (x) < oo, the 
following hold. 


1. The conditional entropy H(xp|x1,...,@n—1) is monotone non-increasing in n. 
2. The order-n entropy is not smaller than the conditional entropy, 
Equation: 


ae) ee Bice s pee 


3. The order-n entropy H,, (z) is monotone non-increasing. 
AL (ay itiig 6 ee) ge od (a ey nny te 


Proof. Part (1): 


Equation: 
bE es se earner Pe SW o Aspe ce eee oe 
= A(e74\ei; soe tnd) 
Part (2): 
Equation: 


[H (a1) + A (xe\a1)+...+H (xp|x1,...,£n—-1)] 


1 
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Part (3): This comes from the fist equality in the proof of (2), because we have the average of a 
monotonely non-increasing sequence. 


Part (4): Both sequences are monotone non-increasing (parts (1) and (3)) and bounded below 
(by zero). Therefore, they both have a limit. Denote H (x) =lim,_,., Hy (x) and 


A (2) lim 40H (a eee): 


Owing to part(2), H> H. Therefore, it suffices to prove H > H, 
Equation: 


1 n+m 


H (ci *) he » H (a;|21,. ‘ rytbe4) 


ps (ones (x) 


A (#j\21, sects 1): 


Now fix n and take the limit for large m. The inequality H < H appears, which proves that 
both limits are equal. 


Coding theorem: ‘Theorem 3 yields for fixed to variable length coding that for a stationary 
source, there exists a lossless code such that the compression rate p,, obeys, 
Equation: 


Ell (x1,...,2n 1 
— [ (21,-.+)%n)] 2H Gh. 
n n 


This can be proved, for example, by choosing / (xz1,...,%n) = [— log P(a1,...,2n)|, which 
is a Shannon code. As 7 is increased, the compression rate p,, converges to the entropy rate. 


We also have a converse theorem for lossless coding of stationary sources. That is, 
Pn = Hy (a) > H. 
Stationary Ergodic Sources 


Consider the sequence x = (...,%_1, £0, £1,...). Let x’ = S, denote a step Vn € Z, 
ZL}, = £41, where S? takes i steps. Let f;, (2) be a function that operates on coordinates 


(Zo,---,@p_1). An ergodic source has the property that empirical averages converge to 
Statistical averages, 
Equation: 
-1 
1 .. @.8.,;N—00 
= SO fk (S;) —— Ef; (2). 
i=0 


In block codes we want 
Equation: 


i; 2 me 


L(@ingis---,€u4n) = H. 


We will be content with convergence in probability, and a.s. convergence is better. 


Theorem 4 Let X be a stationary ergodic source with Hy (x) < oo, then for every 
€ > 0,6 > 0, there exists no (0, €) such that Vn > no (4, €), 
Equation: 


1 _ 
pr {|21(e1,---.0) - H} <e, 
n 


where I (21,...,%n) = — log (Pr (#1,...,2n)). 

The proof of this result is quite lengthy. We discussed it in detail, but skip it here. 
Theorem 4 is called the ergodic theorem of information theory or the ergodic theorem of 
entropy. Shannon (48') proved convergence in probability for stationary ergodic Markov 


sources. McMillan (53') proved L? convergence for stationary ergodic sources. Brieman 
(57'/60') proved convergence with probability 1 for stationary ergodic sources. 


Parametric Models of Information Sources 
In this section, we will discuss several parametric models and see what their entropy rate is. 


Memoryless sources: We have seen for memoryless sources, 
Equation: 


where there are r — 1 parameters in total, 
Equation: 


6 = {p(a),a = 1,2,...,r—1}, 


the parameters are denoted by 6, and a = {1,2,...,r} is the alphabet. 


Markov sources: The distribution of a Markov source is defined as 


Equation: 
uy . 
P (Pi, 895.205 8_) — Plies 4 Le) Il pele) 
i=k+1 
where n > k. We must define {p (a1, @2,---,@k) }(a4,a,...a,)eak initial probabilities and 


transition probabilities, {p (Qpi1 Jat) } There are r* — 1 initial probabilities and (r — 1)r* 


transition probabilities, giving a total of r**+ 


Equation: 


— 1 parameters. Note that 


E {— log p (eile. = At (X;|Xi7;) : ) 8 A(X). 


Therefore, the space of Markov sources covers the stationary ergodic sources in the limit of 
large k. 


Unifilar sources: For unifilar sources, it is possible to reconstruct the set of states that a source 
went through by looking at the output sequence. In the Markov case we have S$; = (X a) , but 


in general it may be more complicated to determine the state. 


To put us on a concrete basis for analysis of unifilar sources, consider a source with M states, 

S = {1,2,..., WM}, and an alphabet a = {1,2,...,r}. In each time step, the source outputs a 
symbol and moves to a new state. Denote the output sequence by x = 21 %2°-- Zn, and the state 
sequence by s = $1S2--- Sn, where s; € S and x; € a. Denote also 

Equation: 


q(s|s') 


Pr{S; = s|St_1 = s'} 
Pr{S; = s|St_4 = s', Si_9, su Y, 


This is a first-order time-homogeneous Markov source. The probability that the next symbol is 
a follows, 
Equation: 


Pr{X; = alS; = s} 
= Pr{X; = al Sy §; X95, 8 215° -}. 


p(als) 


There exists a deterministic function, 
Equation: 


St = g (St-1, Xe-1), 


this is called the next state function. Given that we start at some state S$; = s1, the probability 
for the sequence of states $1, ..., Sn is given by 
Equation: 


p(Xt1S:) = [[ p (80). 


Note the relation 
Equation: 


a(s|s')= >) p(als’). 


a:9(s',a)=s 


To summarize, unifilar sources can be described by a state machine style of diagram as 
illustrated in [link]. 


State machine for selecting the state of a unifilar source. 


Given that an initial state was fixed, a unifilar source with M states and an alphabet of size r 
can be expressed with M(r — 1) parameters. If the initial state is a random variable, then there 
are (MZ — 1 parameters that define probabilities for the initial state, giving 

M(r —1) + _M —1= Mr — 1 parameters in total. In the Markov case, we have M = r*, it is 
a special type of unifilar source. 


Example: 


For the unifilar source that appears in [link], the states can be discerned from the output 
sequence. Let us follow up on this example while discussing more properties of unifilar 
sources. 


Unifilar source for [link]. 


A state S is called transient if once we leave it we cannot go back; Sj is such a state in the 
exaple. A subset & of states is called irreducible if from every state in EF it is possible to move 
to every other state in & in a finite number of steps, and it is impossible to leave EF’. We see that 
{S2, S3} and {.S4} are irreducible sets. The states can be decomposed (in a single unique way) 
into a set of transient states and irreducible components. It can be shown that with probability 1 
the time spent in transient states is finite, and after that the source enters one of the irreducible 
components and stays there forever. Finally, in an ergodic source there may only be one 
irreducible component. Therefre, our example source is not ergodic because it contains two 
such components. 


Given that we are in some state s, which belongs to an irreducible component, the number of 
time steps until we return to s for the first time is a random variable called recurrence time. The 
period of an irreducible component is the largest integer that divides all possible recurrence 
times. 


Example: 
In the unifilar source in [link] the period is 1, because recurrence time could be 2 or 3. 
However, it is impossible for the recurrence time to be 1. 


Unifilar source for [link]. 


An irreducible component with period at least 2 is called periodic. In contrast, a component 
with period 1 is called aperiodic or ergodic, because the probabilities in being in each of the 
states eventually converge to the statistical averages. Note that a periodic component cannot 
yield a stationary distribution. 


An irreducible component F contains asymptotic state probabilities given by the solution of the 


equations, 

Equation: 
Y>a(s'a(s|s’') =4(s), Vs CE. 
siCE 

Additionally, 

Equation: 


lim Pr {S; s|S; sl = qs, (8), 


for aperiodic EZ. The convergence of Pr{S; = s|S; = s'} to q(s) is exponential in t, and 
depends on the eigen-values of the transition probability matrix. In periodic components we do 
not have stationarity, but the following weaker property still holds, 

Equation: 


ee. 
jim eS DPS a8) G55 (8). 


The entropy rate of the output of a unifilar source at state s is given by, 
Equation: 


H (X|S = s) =— 5 p(als)- log p(als). 
acA 
To see this, we begin with a lemma. 


Lemma 1 For a unifilar source, 
Equation: 


H (X,|X?',S1 =s') = S> Pr {S, = |S, = s'}H(X|s). 


ses 


Proof Pr {Xnlxr Si} = Py 1s Sag yi mee Si}, because S,, is deterministic in S,,_1, 
X,,_1. Therefore, owing to the Markov property, 
Equation: 


Pr (X,|X?-*, $1) = Pr(Xn|Sn). 


We can now compute the entropy, 


Equation: 


A(X, | xP SS 3i) — So Pr (XT, Sn)S1) log (Pr (Xn|X7", S1)) 


(XT,Sn) 
= > Pr (xT, Sn| $1) S° EP (Xn|Sn) log (Pr (Xn|Sn)) 
(XT,Sn) Xn 


= » Pr (Sp = s|S1) -H (X|Sp = 8). 


Having proved the lemma, we can move forward with the computation of entropy rate. Note 
that 


Equation: 
1 1i< = 
—H(X, soe ,X,,|S1) Se SoH (Xi|X; : S1) 
n ” =I 
= S°S- Pr (Si = s) H (X|s) 
n l=1 secs 
1 n 
- 3] Fees = 9] mes. 
ses ie I=1 
and 
Equation: 
1 n 
— 5° Pr(S = s) — Foo (s) 
” 1 
Therefore, 
Equation: 


lim —H (X$1S1) = > doo (8)H (X1s). 


N00 
ses 


It is not complicated to show that this is the entropy rate of the unifilar source X, provided that 
it is irreducible and aperiodic. 


Tree Sources: Tree sources are very useful in modeling and processing text and data files in 
general. Consequently, tree sources have been studied in great detail. We will introduce these 
source models here and later in [link] discuss algorithms for universal lossless compression of 
these sources. 


In a tree source, states are arranged as contexts, where a context consists of the last several 
symbols generated by the source. The unique property of a context tree source is that states are 
arranged as leaves of a tree, and the state can be identified by following labels of the branches 
of the tree based on previous symbols, until we reach a leaf. 


Example: 
[link] shows a simple tree source. For this source, a typical input might be 


Equation: 

x = 011001010. 
This input sequence corresponds to the following sequence of states, 
Equation: 


§ = §0950151150505015050150- 


Tree source for [link]. 


As before, for every state (leaf) there exist conditional probabilities for generating the next 
symbol. That said, and maybe surprisingly, a tree source is not necessarily a unifilar source. On 
the one hand, states can be observed from the source symbols. On the other hand, when 
transitioning from a “shallow" leaf to a deep leaf in the tree, the symbol currently being 
generated does not contain enough information to define the new state. To address this problem, 
extra states can be added to the tree. On the other hand, adding extra states to a tree source 
might make the context tree less efficient, because it increases the number of states; see for 
example the detailed discussion of expansion of tree machines by Martin et al. [link]. 


Hidden Markov models: Consider the following state transition probabilities, 
Equation: 


Pr(X; = a, St = s|S4_4 = S's S55 NPG ee ‘) = Pr(X; =a, Sy = s|S_4 = s‘) 
= Pr (a,s|s’). 


The parameter space can be defined, 
Equation: 


0 = {Pr (a, s|s’)}, 


there are M(M - r — 1) parameters in total. 


Note that 
Equation: 
Pr (XT, S?|So) = |] Pr (X:, S:]Sta) 
f=1 
and 
Equation: 


Pr (X?|S0) = 5  Pr(X?, S?|So) 


n 
$1 


S> Il Pr (Xt, S| S¢_-1). 


sy t=1 


It is important to emphasize that the next state is not deterministic. This type of behavior is 
called a hidden Markov model (HMM). 


What are the typical sets for HMM sources? Recall that for an i.i.d. source we saw 
Equation: 


Pr (xq) S|] Bria, 


acA 


this is the main concept underlying the method of types. For a unifilar source, 
Equation: 


Pr (X7) Il Pr (X;|S;) 


Il Pr (als)"*(), 


Ax$S 


For an HMM, there is no such partition into typical sets. 


Universal coding for classes of sources 


We have discussed several parametric sources, and will now start developing 
mathematical tools in order to investigate properties of universal codes that 
offer universal compression w.r.t. a class of parametric sources. 


Preliminaries 


Consider a class A of parametric models, where the parameter set 0 
characterizes the distribution for a specific source within this class, 


{po (-),8 € A}. 


Example: 
Consider the class of memoryless sources over an alphabet a = {1,2,...,r} 
. Here we have 


Equation: 
0 = {p(1), p(2),...,p(r — 1)}. 
The goal is to find a fixed to variable length lossless code that is independent 


of 0, which is unknown, yet achieves 
Equation: 


UXT) m0 


Es Hy (X), 


n 


where expectation is taken w.r.t. the distribution implied by 8. We have seen 
for 
Equation: 


that a code that is good for two sources (distributions) p; and pe exists, 
modulo the one bit loss [link]. As an expansion beyond this idea, consider 
Equation: 


p(#) = | dw (8)po (X), 


where w(6) is a prior. 


Example: 
Let us revisit the memoryless source, choose r = 2, and define the scalar 
parameter 


Equation: 
Then 
Equation: 
po (x) = 0"*) . (1 — ayr*) 
and 
Equation: 


1 
pe= / do "x. (1 — 9), 
0 


Moreover, it can be shown that 
Equation: 


this result appears in Krichevsky and Trofimov [link]. 


Is the source X implied by the distribution p(a) an ergodic source? Consider 
the event lim, +o, - pee. eas oF Owing to symmetry, in the limit of large 
n the probability of this event under p(«) must be =, 

Equation: 


On the other hand, recall that an ergodic source must allocate probability 0 or 
1 to this flavor of event. Therefore, the source implied by p(a) is not ergodic. 


Recall the definitions of pg (x) and p(z) in [link] and [link], respectively. 
Based on these definitions, consider the following, 


Equation: 
Ho (XT) = — >> po(X}) log po (X7) = H(XT|0 = 8), 
XPEA” 
H(Xt) = — > p(X?) log p(X?), 
xX” 


1 


H(X?@) = | dw (0) H(Xp|0 = 8). 


We get the following quantity for mutual information between the random 
variable © and random sequence X i ; 
Equation: 


1(0;X?) = H (Xf) — H(XT|®). 


Note that this quantity represents the gain in bits that the parameter 6 creates; 
more about this quantity will be mentioned later. 


Redundancy 


We now define the conditional redundancy, 
Equation: 


~ [Ey (1(X?)) — Ho (X?)], 


n 


rho = 


this quantifies how far a coding length function / is from the entropy where the 
parameter 6 is known. Note that 
Equation: 


1(Xy') = [ ew (9) Eo (U(Xt)) = A (X74). 


Denote by c, the collection of lossless codes for length-n inputs, and define 
the expected redundancy of a code! € C,, by 
Equation: 
Ry (wl) = | dw(6)rn (1,0), 
A 
R, (w) = inf R, (w,)). 
IEC, 


The asymptotic expected redundancy follows, 
Equation: 


R (w) =lim R;, (w), 


n—- oo 


assuming that the limit exists. 


We can also define the minimum redundancy that incorporates the worst prior 
for parameter, 
Equation: 


R, =sup R,, (w), 
wew 


while keeping the best code. Similarly, 
Equation: 


R’ =lim R,,. 


n—- oo 


Let us derive R,,, 
Equation: 


Ry = supinf [dw (6) — [Bo (I(X1)) ~ H(Xt|© =) 
= supinf —B,|l(X?) — H(X7/®) 


= sup —[H (Xf) — H(Xz|0) 


Ch 
= sup —1(0; Xr) = —-, 


where C’, is the capacity of a channel from the sequence z to the 
parameter [link]. That is, we try to estimate the parameter from the noisy 
channel. 


In an analogous manner, we define 
Equation: 


R; = infsup r, (1,0) 
1 @cA 


po (x”) | 


: 1 
= ae a og es 


1 
= infsup —D(P ; 
nfsup —D(Po||Q) 


where Q is the prior induced by the coding length function l. 


Minimal redundancy 


Note that 
Equation: 
Vw,l, supr,(l,é) > [ 2(a8)r (1,8) 
a A 
& dat / w (d6)rp (1,8). 
lecn A 
Therefore, 
Equation: 


R* =infsup rn (1, 8) >supinf / w (d6)r, (1,0) = R;. 
6 w A 


In fact, Gallager showed that R* = R>. That is, the min-max and max-min 
redundancies are equal. 


Let us revisit the Bernoulli source pg where 6 € A = [0,1]. From the 
definition of [link], which relies on a uniform prior for the sources, i.e., 
w(8) = 1, V6 © A, it can be shown that there there exists a universal code 
with length function / such that 

Equation: 


Nz (1) 


Eo [Il (x")| < nEo hs ( )] + log (n + 1) +2, 


where hz (p) = —p log (p) — (1 — p) log (1 — p) is the binary entropy. That 
is, the redundancy is approximately log (n) bits. Clarke and Barron [link] 
studied the weighting approach, 

Equation: 


p(x) = | dw ()po (2), 


and constructed a prior that achieves R> = R* precisely for memoryless 
sources. 


Theorem 5 [link] For memoryless source with an alphabet of size r, 


6 = (p(0), p(1),---,p(r — 1), 
Equation: 


nk, (w) = Ss : log (=) + | w(as) log (Fer + O, (1), 


where O,, (1) vanishes uniformly as n — oo for any compact subset of A, and 


Equation: 
I(0)2E (- a2 (xi) ) (= = 28 (x;) ) 


is Fisher's information. 


Note that when the parameter is sensitive to change we have large I(@), which 
increases the redundancy. That is, good sensitivity means bad universal 
compression. 


Denote 


Equation: 


vA) 


fi, V8" 


this is known as Jeffrey's prior. Using w(@) = J(@), it can be shown that 
he Sh 


Example: 
Let us derive the Fisher information I(@) for the Bernoulli source, 


Equation: 


po (x) =o" «(1 0)" 
= Inpg(x) =n; (1) Nn d+ n, (0) In (1 - 8) 
O In po (x) 1 1 
apg ag REN) ag 
(- In po (z) ) _ MQ) | m2 (0) _ 2me ()nz (0) 
00 6? (=6)" 6(1 — 6) 


8 In pa (x) \? 1- 2 
oh ( a ‘) | = gt a Og 
1 1 
=@ 19" 
7 1 
~ 6(1—6)- 
Therefore, the Fisher information satisfies I (0) = =-+ 


6(1—6) * 


Example: 
Recall the Krichevsky—Trofimov coding, which was mentioned in [link]. 
Using the definition of Jeffreys’ prior [link], we see that J (0) « Tra 


Taking the integral over Jeffery's prior, 


Equation: 


where we used the gamma function. It can be shown that 
Equation: 


where 
Equation: 
t+1 
pz (x 
Ps (2141/24) = 1 = ) 
ps (24) 
ni; (0) + 5 
Py (41 = O|x;) a a 
nz (1) +3 


Py (41 a 1|z1) = 


Similar to before, this universal code can be implemented sequentially. It is 
due to Krichevsky and Trofimov [link], its redundancy satisfies Theorem 5 by 
Clarke and Barron [link], and it is commonly used in universal lossless 
compression. 


Rissanen's bound 


; eee 4 1 

Let us consider — on an intuitive level — why C,, ~ a aw 

r—1 
2 — 

That is, we would differentiate between each of the r — 1 parameters with ./n 


levels. Now consider a Bernoulli RV with (unknown) parameter 0. 


Expending 


log (7) bits allows to differentiate between (./7) "~" parameter vectors. 


One perspective is that with n drawings of the RV, the standard deviation in 
the number of 1's is O(./7). That is, /7 levels differentiate between 
parameter levels up to a resolution that reflects the randomness of the 
experiment. 


A second perspective is that of coding a sequence of Bernoulli outcomes with 
an imprecise parameter, where it is convenient to think of a universal code in 
terms of first quantizing the parameter and then using that (imprecise) 
parameter to encode the input z. For the Bernoulli example, the maximum 
likelihood parameter 0 yz, satisfies 

Equation: 


Ouy_ =argmax eae! — Jaane 


and plugging this parameter 9 = Oz into pg (x) minimizes the coding length 
among all possible parameters, 8 € A. It is readily seen that 
Equation: 


Nz (1) | 


OuML = 


Suppose, however, that we were to encode with 6’ = Oy, + A. Then the 
coding length would be 


Equation: 


Ig (x) = — log ((6’) a (1 — 6’) a 


It can be shown that this coding length is suboptimal w.r.t. lg,,, (2) by 
n:-O (A?) bits. Keep in mind that doubling the number of parameter levels 


used by our universal encoder requires an extra bit to encode the extra factor 
of 2 in resolution. It makes sense to expend this extra bit only if it buys us at 
least one other bit, meaning that n - O(A?) = 1, which implies that we 


encode 6 yz, to a resolution of 1/4/n, corresponding to O (./n) levels. Again, 
this is a redundancy of + 4 log (n) bits per parameter. 


Having described Rissanen's result intuitively, let us formalize matters. 
Consider {p9,0 € A}, where A C R* is a compact set. Suppose that there 


exists an estimator @ such that 
Equation: 


Yn >n(c): pod | 6(x") — 6 \|> <1 < 5(c), 


where lim. _,.. 6(c) = 0. Then we have the following converse result. 


Theorem 6 (Converse to universal coding [link]) Given a parametric class that 
satisfies the above condition [link], for all « > 0 and all codes / that do not 
know 0, 

Equation: 


K log (n) 


ty (1,0) > (1—«)> os 


except for a class of 6 in B, (n) C A whose Lebesgue volume shrinks to zero 
as m increases. 


That is, a universal code cannot compress at a redundancy substantialy below 
$ log (n) bits per parameter. Rissanen also proved the following achievable 


result in his seminal paper. 


Theorem 7 (Achievable to universal coding [link]) If pg (x) is twice 
differentiable in 0 for every x”, then there exists a universal code such that 


log(n) 
VO EA: rp (1,0) < (1 +e) Ke), 


n 


Universal coding for piecewise i.i.d. sources 


We have emphasized stationary parametric classes, but a parametric class can 
be nonstationary. Let us show how universal coding can be achieved for some 
nonstationary classes of sources by providing an example. Consider 

A= {0,1,...,n} where 

Equation: 


po (x") = Qi (x) - Qe (x51), 


where @ and @» are both know i.i.d. sources. This is a piecewise i.i.d. source; 
in each segment it is i.i.d., and there is an abrupt transition in statistics when 
the first segment ends and the second begins. 


Here are two approaches to coding this source. 


1. Encode the best index @yyz, using [log (n + 1)] bits, then encode 
P61, (£”). This is known as two-part code or plug-in; after encoding the 
index, we plug the best parameter into the distribution. Clearly, 
Equation: 
(x) = min [— log po (x)| + [log (n + 1)] 


0<0<n 


< — log po (x)+ log (n +1) +2. 


2. The second approach is a mixture, we allocate weights for all possible 
parameters, 
Equation: 


We) = — be (2 Sn) 


1 
— log (Pin ")) 


— log (Dor (x))+ log (n aE 1). 


x 


Merhav [link] provided redundancy theorems for this class of sources. 
Algorithmic approaches to the mixture appear in Shamir and Merhav [link] 
and Willems [link]. 


The theme that is common to both approaches, the plug-in and the mixture, is 
that they lose approximately log (n) bits in encoding the location of the 
transition. Indeed, Merhav showed that the penalty for each transition in 
universal coding is approximately log (7) bits [link]. Intuitively, the reason 
that the redundancy required to encode the location of the transition is larger 
than the + log (n) from Rissanen [link] is because the location of the 
transition must be described precisely to prevent paying a big coding length 
penalty in encoding segments using the wrong i.i.d. statistics. In contrast, in 
encoding our Bernoulli example an imprecision of = in encoding @y¢z in the 


Vn 
first part of the code yields only an O(1) bit penalty in the second part of the 


code. 


It is well known that mixtures out-compress the plug-in. However, in many 
cases they do so by only a small amount per parameter. For example, Baron et 
al. showed that the plug-in for i.i.d. sources loses approximately 1 bit per 
parameter w.r.t. the mixture. 


Universal compression for context tree sources 


Semi-predictive approach 


Recall that a context tree source is similar to a Markov source, where the 
number of states is greatly reduced. Let T be the set of leaves of a context 
tree source, then the redundancy is 


Equation: 
r< Air) (108 (Fz) +o), 


where |7' is the number of leaves, and we have log (4) instead of 


oa 


| 
the redundancy for a Markov representation of the tree source 7’ is much 
larger. Therefore, tree sources are greatly preferable in practice, they offer a 
significant reduction in redundancy. 


log (n), because each state generated symbols, on average. In contrast, 


How can we compress universally over the parametric class of tree sources? 
Suppose that we knew T,, that is we knew the set of leaves. Then we could 
process x sequentially, where for each x; we can determine what state its 
context is in, that is the unique suffix of i that belongs to the set of leaf 
labels in 7’. Having determined that we are in some state s, 

Pr (x; = O|s, ai") can be computed by examining all previous times that 
we were in state s and computing the probability with the Krichevsky- 
Trofimov approach based on the number of times that the following symbol 
(after s) was 0 or 1. In fact, we can store symbol counts n, (s,0) and 

Nz (s,1) for all s € T, update them sequentially as we process x, and 
compute Pr (a = O|s, a.) efficiently. (The actual translation to bits is 
performed with an arithmetic encoder.) 


While promising, this approach above requires to know 7’. How do we 
compute the optimal T’” from the data? 


Tes ge Ee 


Tree pruning in the semi-predictive 
approach. 


Semi-predictive coding: The semi-predictive approach to encoding for 
context tree sources [link] is to scan the data twice, where in the first scan 
we estimate J’ and in the second scan we encode x from T' . as described 
above. Let us describe a procedure for computing the optimal ‘i among 
tree sources whose depth is bounded by D. This procedure is visualized in 
[link]. As suggested above, we count n,, (s, a), the number of times that 
each possible symbol appeared in context s, for all s € a”, a € a. Having 
computed all the symbol counts, we process the depth-D tree in a bottom- 
top fashion, from the leaves to the root, where for each internal node s of 
the tree (that is, s € a? where d < D), we track 7 the optimal tree 
structure rooted at s to encode symbols whose context ends with s, and 
MDL(s) the minimum description lengths (MDL) required for encoding 
these symbols. 


Without loss of generality, consider the simple case of a binary alphabet 
a = {0,1}. When processing s we have already computed the symbol 
counts nz (0s, 0) and nz (0s, 1), nz (15,0), nz (1s, 1), the optimal trees 


ee and T.; and the minimum description lengths (MDL) MDL(0s) and 
MDL(1S). We have two options for state s. 


1, KeepT yg and ioe The coding length required to do so is 
MDL(0S) + MDL(1S) + 1, where the extra bit is spent to describe 
the structure of the maximizing tree. 

2. Merge both states (this is also called tree pruning). The symbol counts 
will be nz (8, a) = nz (08, a) + nz (1s, a),a € {0,1}, and the 
coding length will be 
Equation: 


KT(nz (s,0), nz (s,1)) +1, 


where K'T(-, -) is the Krichevsky-Trofimov length [link], and we again 
included an extra bit for the structure of the tree. 


We note in past that there is no need to spend a bit to encode leaves of depth 
D. To see this, consider a procedure for encoding the structure of a tree: 


Example: 

Consider the tree sourced depicted in [link]. In order to encode the 
structure of this tree, we will utilize the following procedure. (Such a 
procedure has appeared, for example, in [link].) 


Tree used in [link] to demonstrate how the structure 
of the source is encoded. 


Start from root. [procedure(root)] 
1. If node S is of depth D (maximum), then return. 
2. If node S is internal node, then { 
encode 0 
procedure(0S) 
procedure(1S) 
} else encode 1. 
3. return. 


Let us now simulate the procedure, the procedure will traverse through the 
following states of the tree in [link] while outputting the corresponding bits. 


Source root 0 1 01 O01 101 11 


Encoded 
symbol 


Returning to tree pruning, following [link] we see that we must initialize 
MDL (s) = KT(n, (s,0), nz (s, 1)) for s of full depth |s| = D without 
the extra bit. 


At the end of the pruning procedure, T the maximizing tree for the root, 
will be the optimal tree for universal coding. 


Burrows wheeler Transform 


The Burrows Wheeler transform (BWT) was proposed by Burrows and 
Wheeler in 1994 [link] (see also the analysis by Effros et al. [link] and 
references therein). It is an invertible permutation sort that sorts symbols 
according to their contexts. That way, the symbols that were generated by 
the same state of the context tree are grouped together, which as we will see 
is advantageous. 


To compute the BWT, we first compute all cyclical shifts of the input x. 
Next, we sort the cyclical shifts. The output of the BWT consists of y, the 
last column of the matrix of sorted shifts, and 2 the index of the original 
version. We illustrate with an example. 


Example: 
Consider the input x = banana. First, we compute the cyclic shifts and 
their sorts. 


All Shifts Sorted 


banana abanan 
abanan anaban 
nabana ananab 
anaban banana 
nanaba nabana 
ananab nanaba 


The output of the BWT consists of y = nnbaaa, the last column of the 
matrix of sorted shifts (to the right), and the index 7 = 4 containing the 
original input. 


Interestingly, we can recover x from y and 2. Seeing that y is structured and 


thus quite compressible, the BWT can be used as a compression system; a 
building block that illustrates such a system appears in [link]. 


Y 


Typical compression system using the Burrows Wheeler 
transform [link]. 


To see that the BWT is invertible, let us work out how to do this by 
continuing our example. 


Example: 
In the matrix of sorted shifts, column 1 is a sorted version of column n, 
which we know. 


Column 1 Column n 
a n 
a n 
a b 
b a 
n a 
n a 


Now take column n and put it before column 1: 


Column n Column 1 


n a 
b a 
a b 
a n 
a n 


We now sort these rows, which each consist of 2 symbols: ab, an, an, ba, 
na, and na. Now fill column 2 of the sorted shifts matrix accordingly. 


Columns 1-2 Column n 
ab n 
an n 
an b 
ba a 
na a 
na a 


The entire matrix can be unraveled, and the row containing the original x is 
indexed by 7. 


What is the BWT good for? The key property of the BWT is that symbols 
generated by the same state are grouped together in y. To see this, note how 
the last column n can be rotated to a position to the left of column 1, and 
symbols that came before the same prefix appear together. (To bunch 
together symbols generated by the same suffix, we can reverse the order of 
symbols in « before running the BWT.) Therefore, y has the form of a 
piecewise i.i.d. sequence [link], where segments generated by the same 
State of the context tree are bunched together. 


As a consequence of the structure of y, it is easy to see that it can be 
compressed with the following redundancy, 
Equation: 


r< cai (108 (a + 0(1)) + |T| log (n), 


where the new |7'| log (n) term arises from coding the locations of 
transitions between segments (states of the tree) in the BWT output. Not 
only is the BWT convenient for compression, but it is amenable to fast 
computation. Both the BWT and its inverse can be implemented in O(n) 
time. This combination of great compression and speed has made the BWT 
quite popular in compressors that have appeared since the late 1990s. For 
example, the bzip2 archiving package is very popular among network 
administrators. 


That said, from a theoretical perspective the BWT suffers from an 
extraneous redundancy of |T'| log (n) bits. Until this gap was resolved, the 
theoretical community still preferred the semi-predictive method or another 
approach based on mixtures. 


Semi-predictive coding using the BWT 


Another approach for using the BWT is to use y only for learning the MDL 
tree source T’. To do so, note that when the BWT is run, it is possible to 
track the correspondences between contexts and segments of the BWT 
output. Therefore, information about per-segment symbol count is 


available, and can be easily applied to perform the tree pruning procedure 
that we have seen. Not only that, but some BWT computation algorithms 
(e.g., suffix tree approaches) maintain this information for all context 
depths and not just bounded D. In short, the BWT allows to compute the 
minimizing tree T’ in linear time [link]. 


Given the minimizing tree T ’, it is not obvious how to determine which 
state generated each character of y (respectively, x) in linear time. It has 
been shown by Martin et al. [link] that this step can also be performed in 
linear time by developing a state machine whose states include the leaves of 
T. The result is a two part code, where the first part computes the optimal 
T° via BWT, and the second part actually compresses x by tracking which 
state of T" generated each of the symbols. To summarize, we have a linear 
complexity algorithm for compressing and decompressing a source while 
achieving the redundancy bounds for the class of tree sources. 


Context Tree Weighting 


We discussed in [link] for the problem of encoding a transition between two 
known i.i.d. distributions that 
Equation: 


Pala ) > = max {po, (2)}, 


Therefore, a mixture over all parameter values yields a greater probability 
(and thus lower coding length) than the maximizing approach. Keep in 
mind that finding the optimal MDL tree source T is analogous to the plug- 
in approach, and it would reduce the coding length if we could assign the 
probability as a mixture over all possible trees, where we assign trees with 
fewer leaves a greater weight. That is, ideally we want to implement 
Equation: 


Pr (2) = J) 2. pe (a), 


T 


where |code(T’)| is the length of the encoding procedure that we discussed 
for the tree structure JT’, and pr () is the probability for the sequence x 
under the model 7’. 


Willems et al. showed how to implement such a mixture in a simple way 
over the class of tree sources of bounded depth D. As before, the algorithm 
proceeds in a bottom up manner from leaves toward the root. At leaves, the 
probability p, assigned to symbols that were generated within that context s 
is the Krichevsky-Trofimov probability, px (s, x) [link]. For s that is an 
internal node whose depth is less than D, the approach by Willems et 

al. [link] is to mix (i) the probabilities of keeping the branches for Os and 1s 
and (ii) pruning, 

Equation: 


1 1 
Ps = gPKh (s, x) “1 9 Pos *Pis- 


It can be shown that this simple formula allows to implement a mixture 
over the class of bounded depth context tree sources, thus reducing the 
coding length w.r.t. the semi-predictive approach. 


In fact, Willems later showed how to extend the context tree weighting 
(CTW) approach to tree sources of unbounded depth [link]. Unfortunately, 
while the basic bounded depth CTW has complexity that is comparable to 
the BWT, the unbounded CTW has potentially higher complexity. 


Beyond lossless compression 


Universal lossy compression 


Consider x € a”. The goal of lossy compression [link] is to describe %, also of length n but 
possibly defined over another reconstruction alphabet @ 4 a, such that the description requires 
few bits and the distortion 

Equation: 


is small, where d(-, -) is some distortion metric. It is well known that for every d(-, -) and 
distortion level D there is a minimum rate R(D), such that & can be described at rate R(D). The 
rate R(D) is known as the rate distortion (RD) function, it is the fundamental information 
theoretic limit of lossy compression [link], [link]. 


The invention of lossy compression algorithms has been a challenging problem for decades. 
Despite numerous applications such as image compression [link], [link], video compression [link], 
and speech coding [link], [link], [link], there is a significant gap between theory and practice, and 
these practical lossy compressors do not achieve the RD function. On the other hand, theoretical 
constructions that achieve the RD function are impractical. 


A promising recent algorithm by Jalali and Weissman [link] is universal in the limit of infinite 
runtime. Its RD performance is reasonable even with modest runtime. The main idea is that the 
distortion version % of the input x can be computed as follows, 

Equation: 


© = argmin {He (w") — Bd (2",w")}, 


wrean 


where 8 < 0 is the slope at the particular point of interest in the RD function, and Hj, (w”) is the 
empirical conditional entropy of order k, 
Equation: 


thy (ua 
Hy, (w") = = dm (u*, a) log (<=“3,,). 


a,u alca Tw (ak, a’) 


* is a context of order k, and as before ny, (u*, a) is the number of times that the symbol 


a appears following a context u* in w”. Jalali and Weissman proved [link] that when 
k = o(log (n)), the RD pair (4 (a7 \0 (2", a”)) converges to the RD function 


asymptotically in n. Therefore, an excellent lossy compression technique is to compute % and then 
compress it. Moreover, this compression can be universal. In particular, the choice of context 


where u 


order k = o(log (m)) ensures that universal compressors for context tress sources can emulate the 


coding length of the empirical conditional entropy H k (a), 


Despite this excellent potential performance, there is still a tremendous challenge. Brute force 


computation of the globally minimum energy solution z” involves an exhaustive search over 
exponentially many sequences and is thus infeasible. Therefore, Jalali and Weissman rely on 
Markov chain Monte Carlo (MCMC) [link], which is a stochastic relaxation approach to 
optimization. The crux of the matter is to define an energy function, 

Equation: 


e(w”) = Hy (w”) — Bd(a",w"). 


The Boltzmann probability mass function (pmf) is 
Equation: 


where t > 0 is related to temperature in simulated annealing, and Z; is the normalization 
constant, which does not need to be computed. 


Because it is difficult to sample from the Boltzmann pmf [link] directly, we instead use a Gibbs 
sampler, which computes the marginal distributions at all n locations conditioned on the rest of 
w” being kept fixed. For each location, the Gibbs sampler resamples from the distribution of w; 
conditioned on w”\' & {w,, : 2 #1} as induced by the joint pmf in [link], which is computed as 
follows, 
Equation: 
Pr (w; = b, w"\’) 

Pr (w"\*) 

Pr (w; = b, w"\') 
Nivea Pr (wi = ',w™‘) 

z exp {—te (w; = b, w"\") 

ye Z exp {+e (we bia) | 


Pr (wi = b|w") = 


We can now write, 
Equation: 


e(yi=b,y") =e (yi =a) + AH (6, a) ~ BAd(b, a), 


where AH; (y' "by. 4, a) is the change in Hj, (w”) when y; = a is replaced by b, and 
Ad (b, a, x;) = d(b, x;) — d(a,;) = (b— 2)? — (a — 2)? is the change in distortion. 
Combining [link] and [link], 

Equation: 


r (wi =blw"’) = 
3 ( b ) Sy exp {—F[e (a) + AH (b',a) — BAd(b’,a)]} 


exp {7 0} 


Sy exp {-F[AH (b',a) — BAd(b/,a) — AH (b, a) + BAd (b, a)]} 
1 


Dy exp {—1[AH (b',b) — BAd(b’,d)]} 


The maximum change in the energy within an iteration of MCMC algorithm is then bounded by 
Equation: 


A; =max max max In AH; (0, b’) + c4Ad (6, b’) |. 


1<i<nweanrd,b'ea 


We refer to the resampling from a single location as an iteration, and group the n possible 
locations into super-iterations.[footnote] 

Baron and Weissman [link] recommend an ordering where each super-iteration scans a 
permutation of all n locations of the input, because in this manner each location is scanned fairly 
often. Other orderings are possible, including a completely random order. 


During the simulated annealing, the temperature ¢ is gradually increased, where in super-iteration 
i we use t = O(1/ log (2)) [link], [link]. In each iteration, the Gibbs sampler modifies w” in a 
random manner that resembles heat bath concepts in statistical physics. Although MCMC could 
sink into a local minimum, we decrease the temperature slowly enough that the randomness of 
Gibbs sampling eventually drives MCMC out of the local minimum toward the set of minimal 
energy solutions, which includes z”, because low temperature t favors low-energy w”. Pseudo- 
code for our encoder appears in Algorithm 1 below. 


Algorithm 1: Lossy encoder with fixed reproduction alphabetInput:2” € a”, a, 8, c, rOutput: bit- 
stream Procedure: 


1. Initialize w by quantizing x with @ 
2. Initialize ny (-, -) using w 
3. for r = 1 to R do // super-iteration 


cnA. 
tp — fost) // temperature 


Draw permutation of numbers {1, ...,n} at random 
for r' = 1tondo 

Let 7 be component t’ in permutation 

Generate new w; using f, (w; = -|w"\’) 


ONaw - 


9, Update ny (-, -) [+] 
10. Apply CTW to w” // compress outcome 


Computational issues 
Looking at the pseudo-code, it is clear that the following could be computational bottlenecks: 


1. Initializing n (-,-) - a naive implementation needs to scan the sequence w” (complexity 
O(n)) and initialize a data structure with |@ ine elements. Unless K Slog) (n), this is 
super linear in n. Therefore, we recall that K = o(log (n)), and initializing n,, (-, -) requires 
linear complexity O(n). 

2. The inner loop is run rn times, and each time computing Pr (w j= blw ) for all possible 
b € a might be challenging. In particular, let us consider computation of Ad and AH. 


1. Computation of Ad requires constant time, and is not burdensome. 

2. Computation of AH requires to modify the symbol counts for each context that was 
modified. A key contribution by Jalali and Weissman was to recognize that the array of 
symbol counts, ny (-,-), would change in O(k) locations, where k is the context order. 
Therefore, each computation of AH requires O(k) time. Seeing that || such 
computations per iteration are needed, and there are rn iterations, this is O(krn|@|). 

3. Updating n,, (-,-) after w; is re-sampled from the Boltzmann distribution also requires 
O(k) time. However, this step is performed only once per iteration, and not ja | times. 
Therefore, this step requires less computation than step (b). 


Lossy compression of an analog input 


So far, we have described the approach by Jalali and Weissman to compress x” € a”. Suppose 
instead that x € IR” is analog. Seeing that MCMC optimizes w” over a discrete alphabet, it would 
be convenient to keep doing so. However, because @ is finite, and assuming that square distortion 
is used, i.i., d (xi, wi) = (ai — wi)’, we see that the distortion could be large. 


Fortunately, Baron and Weissman showed [link] that the following quantizer has favorable 
properties, 
Equation: 


2 2 2 

x ene ees eae 1 ty 
oe en ’ 
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where 7¥ is an integer that grows with n. That is, as yy grows the set @ of reproduction levels 
quantizers a wider internal more finely. Therefore, we have that the probability that some z; is an 
outlier, ie., |2;|< -y, vanishes with n, and the effect of outliers on d (x, Z) vanishes. Moreover, 
because the internal is sampled more finely as y increases with n, this quantizer can emulate any 
codebook with continuous levels, and so in the limit its RD performance converges to the RD 
function. 


Rate distortion performance 


To evaluate the RD performance, we must first define the RD function. Consider a source X that 
generates a sequence z” € R”. A lossy encoder is a mapping e : R” + @, where @ C R” isa 
codebook that contains a set of codewords in R”. Let @” (x, D) be the smallest cardinality 
codebook for inputs of length n generated by X such that the expected per-symbol distortion 
between the input 2” and the codeword e (x) € @” (x, D) is at most D. The rate R” (x, D) is 
the per-symbol log-cardinality of the codebook, 

Equation: 


R (e, D) = — log (|€" (0, D)|), 


This is an operational definition of the RD function in terms of the best block code. In contrast, it 
can be shown that the RD function is equal to a minimization over mutual information, yielding a 
different flavor of definition [link], [link]. 


Having defined the RD function, we can describe the RD performance of the compression 
algorithm by Baron and Weissman for continuous valued inputs. 


Theorem 8 Consider square error distortion, d (wi, Zi) — (x; — @;) . let X be a finite variance 
stationary and ergodic source with RD function R(X, D), and use the MCMC algorithm with data 
independent reproduction alphabet @ and temperature that decays sufficiently slowly. Let w?’ be 
the approximation to x” after r super-iterations. Then the length of context tree weighting (CTW) 
applied to w?’ converges as follows, 

Equation: 


1 7 
lim lim E|—|CTW (W,")| — 6d (2",W2 in [R(X,D) — BD). 
im lim FE) — | (W,")| — Bd (2", W,"))n — 00 min [R(X, D) — BD] 


n—-oor—-oo D>0 


Remark 1 Let us make some comments. 


¢ This result implies that CTW is used to compress w;,. after it has been generated by MCMC. 
Instead of CTW, other universal compression approaches can be used. 

e The same result can be proved for /,, distortion metrics, and not just lo. 

e A result with similar flavor was proved by Jalali and Weissman for discrete alphabets [link]. 


Adaptive reproduction levels: The above algorithm is promising from a theoretical perspective, 
but is of limited practical interest. In order to approach the RD function closely, @ may need to be 
large, which slows down the algorithm. 


Our approach to overcome the large alphabet is inspired by the work by Rose [link], who showed 
that in many cases the optimal RD codebook uses a small discrete-valued reproduction alphabet. 
This runs counter to our intuition that a continuous reproduction alphabet is needed to compress a 
continuous valued source. Therefore, we propose an algorithm that allows for reduction of the 
alphabet size while supporting adaptive reproduction levels. 


We map the input x” to a sequence z” over a finite alphabet, where actual numerical values are 
obtained by a scalar function, g(z;). Ideally, g(-) should minimize distortion. Because we focus on 
square error distortion, the optimal g* (-) is conditional expectation, 

Equation: 


yess vi 
a 1 


g (a) = E[z,|z; =a] = 


Keep in mind that in an actual compression system, the encoder knows x”, but the decoder does 
not. Therefore, the encoder must describe (a quantized version of) g* (-) to the decoder. Putting 
aside these details, the main point is that the adaptive reproduction alphabet algorithm also 
acheives the RD function for stationary and ergodic sources while allowing a reduction in the size 
of the alphabet, thus accelerating the algorithm. 


Universal signal reconstruction 


Scalar channel: Consider a scalar channel, y = x + z, where x, y, z € IR”, x is generated by a 
stationalry and ergodic source X, and z; ~ N (0,1) is zero mean unit norm Gaussian noise. To 
avoid taking the analysis toward continuous valued x, we can simplify the setting and temporarily 
consider discrete-valued x, y, and z; for example, we can consider a quantized version of the 
continuous problem. 


Donoho [link] proposed the following reconstruction, 
Equation: 


= argmin K(w), 
ws.t.||Y—w||3<n 


where w is a possible reconstruction sequence, and K(w) is the Kolmogorov complexity of w. 
Donoho calls this reconstruction algorithm the Kolmogorov sampler. 


What is the Kolmogorov complexity [link], [link], [link] of the sequence w? It is the length of the 
shortest computer program that outputs w and then halts. One may ask what computer language 
we can use for the program; after all, it is easy to imagine some programming language being 
well-tuned to a specific sequence w. However, concepts such as universal Turing machines [link] 
can be used to put forward an abstract programming language that will run on any hypothetical 
computer. The key point is that the computer is limited to a state machine that processes an 
infinite-length tape, and one state machine can emulate another by programming a compiler onto 
the tape. 


Is it practical to talk about the shortest program that outputs w and then halts? Not really, because 
we would need to run the Turing machine [link] (the computer) on an exponential number of 
programs. Moreover, it is well known that for some programs the Turing machine never halts 
(there are endless loops), and technical conditions such as time-out mechanisms will need to be 
added in. Therefore, computing the Kolmogorov complexity A (w) is impractical [link], [link], 


[link]. However, for w that was generated by a stationary ergodic source, it is possible to 
approximate K'(w) by H;, (w), the empirical entropy [link], [link]. 


What is the performance of Donoho's Kolmogorov sampler? Donoho proved that % has the 
following distribution, X ~ p (X|Y) [link]. That is, % is sampled from the posterior. Therefore, 
given y, we have that & and z are both generated by the same distribution. 


Consider now that %,y;7spz, the minimum mean square error estimator, is the center of mass of the 
posterior. Seeing that  — Zymse and x — Xywuse are orthogonal, 
Equation: 


E||| « — 2 ||3| = Ell] «Fase lo] + Ell @umse — = ||]; 


and so E||| « —  ||3], the mean square error, is double the MMSE [link]. 


Let us pause to reflect about this result. When the SNR is high, ie., || « ||} >> || z ||3, the MMSE 
should be rather low, and double the MMSE seems pretty good. On the other hand, when the SNR 
is low, the MMSE could be almost as large as E||| « ||3], and double the MMSE could be larger — 
as much as twice larger — than E[|| z ||3]. That is, gussing @ = E [2] could give better signal 
estimation performance than using the Kolmogorov sampler. This pessimistic result encourages us 
to search for better signal reconstruction methods. 


Arbitrary channels: So far we considered the Kolmogorov sampler for the white scalar channel, 
y = x + z. Suppose instead that x is processed or measured by a more complicated system, 
Equation: 


y= J(x) +z. 


Note that J is known, e.g., in a compressed sensing application [link], [link] would be a known 
matrix. An even more involved system would be y = (J @ x) © z, where J © z is application of 
a mapping to x, J @ x € R®, and @z denotes application of a random noise operator to J ® x. To 
keep the presentation simple, we use the additive noise setting [link]. 


How can the Kolmogorov sampler [link] be applied to the additive noise setting? Recall that for 
the scalar channel, the Kolmogorov sampler minimizes K(w) subject to || y — w ||} <n. For the 
arbitrary mapping J with additive noise [link], this implies || y — J (w) ||} <n. Therefore, we 


get 
Equation: 


f= argmin K(w). 


s.t.|Y¥—J(w)||3<n 


Another similar approach relies on optimization via Lagrange multipliers, 
Equation: 


© =argmin K (w)— log, (f.(y— J (w))), 


where the lagrange multiplier is 1, because both K(w) and log, (fz (y— J (w))) are quantified 
in bits. 


What is the performance of the Kolmogorov sampler for an arbitrary J? We speculate [link], 
[link] that Z is generated by the posterior, and so F | xr ia] is double the MMSE, where 


expectation is taken over the source X and noise z. These results remain to be shown rigorously. 


Convergence of MCMC algorithm 


We will now prove a substantial result — that the MCMC algorithm [link], [link], [link] converges 
to the globally minimal energy solution for the specific case of compressed sensing [link], [link]. 
An extension of this proof to arbitrary channels J is in progress. 


R™*n 


If the operator J in [link] is a matrix, and we denote it by ® € | where m < n, then the 
setup is known as compressed sensing (CS) [Link], [link] and the estimation problem is commonly 
referred to as recovery or reconstruction. By posing a sparsity or compressibility requirement on 
the signal and using it as a prior during recovery, it is indeed possible to accurately estimate x 
from y in the CS setting. 


With the quantization alphabet definition in [link], @ will quantize x to a greater resolution as N 
increases. We will show that under suitable conditions on fx, performing maximum a posteriori 
(MAP) estimation over the discrete alphabet @ asymptotically converges to the MAP estimate 
over the continuous distribution fx. This reduces the complexity of the estimation problem from 
continuous to discrete. 


We assume for exposition that we know the input statistics f. Given the measurements y, the 
MAP estimator for x has the form 
Equation: 


MAP = argmax fx (w) fyix (ylw). 


Because z is i.i.d. Gaussian with mean zero and known variance a7, 
Equation: 


fy\x (y|w) = eye cally-Pu” 


My 1 
207, 
Plugging into [link] and taking log likelihoods, we obtain 
Equation: 


_—M/2 , 
where cy = (270%) and cy = are constants and || - || denotes the Euclidean norm. 


tap = argminY~ (w), 
WwW 


where &~ (-) denotes the objective function (risk) 
Equation: 


U* (w) = — In (fx (w)) + e2|| y — Sw ||’; 


our ideal risk would be ©* (x yy4p). 


Instead of performing continuous-valued MAP estimation, we optimize for the MAP in the 
discretized domain @. We begin with a technical condition on the input. 


Condition 1 [link] We require that the probability density has bounded derivatives 
Equation: 


d 


Ln, 


In (fx (x))| <p, 


where _f is the derivative w.r.t. entry n of z,n € {1,...,N}, and p > 0. 


Ln, 


Let %,zap be the quantization bin in @ nearest to x 4p. Condition 1 ensures that a small 
q Condition 1 


perturbation from xy 4p to Zap does not change In (fx (-)) by much. We use this fact to prove 
that b* (Zap) is sufficiently close to y* (xmap) asymptotically. 


Theorem 9 [link] Let 6 € R“** be ani.id. Gaussian measurement matrix where each entry has 
mean zero and variance ae Suppose that Condition 1 holds and the aspect ratio d = “ > 0, and 
let the noise z € R™ be iid. zero-mean Gaussian. Then for all € > 0, the quantized estimator 
Zu AP Satisfies 

Equation: 


yx (tap) < yx (Zap) < Ut (aap) + Ne 


almost surely as N — oo. 

Given that Theorem 9 shows that the risk penalty due to quantization vanishes asymptotically in n 
, we now describe a Kolmogorov-inspired estimator for CS over a quantized grid. We define the 
energy e(w”): 

Equation: 


e(w") & nH (w") + call y™ — Su" ||’, 


where C4 = C2 log, (e). 


Ideally, our goal is to compute the globally minimum energy solution 
Equation: 


ae = argmine (w”). 


wrean 


Theorem 10 [link] Let X be a bounded stationary ergodic source. Then the outcome wy?’ of 
Algorithm 1 after r iterations obeys 
Equation: 


H 
Jim ¢ (wr) a oy la eam) 


We define a stochastic transition matrix P(,) from @” to itself given by the Boltzmann distribution 
for super-iteration r. Similarly, 7r(,) defines the stable state distribution on a” for P(r), satisfying 
T(r) P(r) = Mr): 


Definition 4 [link]The Dobrushin's ergodic coefficient of a Markov chain transition matrix P is 
denoted by 6(P) and defined as 
Equation: 


6(P) = max ea 


4 max ll pi — 25 ll 


where p; denotes the i” row of P, l1<n< N. 


From the definition, 0 < 6(P) < 1. Moreover, the ergodic coefficient can be rewritten as 
Equation: 
N 


6(P)=1- min ae (Dik Pik) 


where p;; denotes the entry of P at the i*” row and j column. 


We group the product of transition matrices across super-iterations as 
Equation: 


Porsr2) = II Poy: 


T=T1 


There are two common characterizations for the stable-state behavior of a non-homogeneous MC. 


Definition 5 [link] A non-homogeneous MC is called weakly ergodic if for any distributions ps 
and v over the state space Y, and any r; € N, 
Equation: 


limsup,, -5.0| LP 1-999) _ VP m1 92) ll, = 0. 


Similarly, a non-homogeneous MC is called strongly ergodic if there exists a distribution 7 over 
the state space Y such that for any distribution yz over , and any r; € N, 
Equation: 


limsup,,, 5.0 || LP(r, 19) _T ll, = 0. 


We will use the following two theorems from [link] in our proof. 


Theorem 11 [link] A Markov chain is weakly ergodic if and only if there exists a sequence of 

integers 0 < ry < rq <... such that 

Equation: 
Co 
>» 0 = 5(Porisriss))) ee 
i=1 
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Theorem 12 [link] Let a Markov chain be weakly ergodic. Assume that there exists a sequence of 
probability distributions {a(n ben on the state space .Y such that 7(,) P(r) = m(r). Then the 
Markov chain is strongly ergodic if 

Equation: 


CO 


ll T(r) — Tera) I], < c- 
r=1 


The rest of proof is structured as follows. First, we show that the sequence of stable state 
distributions for the Markov chains (MC) used by the MCMC algorithm converges to a uniform 
distribution over the set of sequences that minimize the energy function as the iteration count t 
increases. Then, we show using Theorem 11 and Theorem 12 that the non-homogeneous MC used 
in the MCMC algorithm is strongly ergodic, which by the definition of strong ergodicity implies 
that MCMC always converges to the stable distribution found above. This implies that the 
outcome of the MCMC algorithm converges to a minimum-energy solution as t —> oo, 
completing the proof of Theorem Theorem 10. 


We therefore begin by finding the stable state distribution for the non-homogeneous MC used by 
the MCMC algorithm. At each super-iteration 7, the distribution defined as 
Equation: 


exp (- re (w")) 


Yo amcan EXP (- ee (2")) 
1 


oancan xp (—E(€(2") — e(w"))) 


I 


T(r) (w”) 


satisfies 7 (p)P(p) = 7p). We can show that the distribution 7, converges to a uniform 
distribution over the set of sequences that minimize the energy function, i.e., 
Equation: 


, 0 w¢é ZH, 
Ammo =) 1 um x0, 


where # = {w" € @" s.t. €(w”) =min»egn €(z") }. To show [link], we will show that 
T(r) (w”) is increasing for w” € # and eventually decreasing for w” € # C. Since for 
w" € # and z” € a”, e(z”) — e(w”) > 0, and so for r1 < rz we have 

Equation: 


Y exp (-F-€e(@")-e(w"))) > exp (- Fe") -ew"y), 


2zEQ” wea” 


which together with [link] implies 7(,,) (ao) =< T (r3) (w”). On the other hand, if w” ¢ #, then 
Equation: 
1 
T(r) (w”) 


= Y ew{-Zee-ew'y} 


Z”:e(z”) >e(w”) 


ZMe(z™ 


For sufficiently small t, the right hand side (more precisely, the second and third lines) of [link] is 
dominated by the second term (third line), which increases when ¢ decreases, and therefore 
T(r) (w™) decreases for w” € 7 as t increases. Finally, since all sequences w” € # have the 


same energy €(w”), it follows that the distribution is uniform over the symbols in #. 


Having shown convergence of the non-homogenous Markov chain's stable state distributions, we 
now show that the non-homogeneous MC is strongly ergodic. The transition matrix P/,) of the 
MC at iteration t depends on the temperature ¢ used within MCMC algorithm. We first show that 


the MC used in the MCMC algorithm is weakly ergodic via Theorem 11; the proof of the 
following Lemma is given at the end of this section. 


Lemma 2 The ergodic coefficient of P(,) for any r > 0 is upper bounded by 
Equation: 


(Pry) < 1-— exp {onde}, 


where A is defined in [link]. 


Let w}, ws be two arbitrary sequences in a”. The probability of transitioning from a given state 
to a neighboring state in an iteration within iteration r’ of super iteration r of the MCMC 
algorithm is given by [link], and can be rewritten as 

Equation: 


ce (wr ~ b[w"\”’) = EE (w = b|w"\” ) 


1 r'—1 n r'—1 n 
exp (-(¢ (wi bur) = Euan (w w41))) 
i 
t 


Vive exp (— 7 (€ (wh "bw 1) — émings (WE, wF41))) 


where €min,r’ (wi ~* bw” al =MiNy cg € (wi bw”, s) Therefore the smallest probability of 


transition from w? to w7 within super-iteration r is bounded by 
Equation: 


wy wyean 


min P(r) (wz|w7) > TI al = ay" 


Using the alternative definition of the ergodic coefficient given in [link], 
Equation: 


(P(r) =1-— min S° min (Por) (z"|w7), Pr) (z"|w3)) 


w" wrean = 
172 2ZrEqn 
n 


exp (—inAr) 1 
1—|a@ a = 1- exp (~j-nar). 


[a 


Using Lemma 2, we can evaluate the sum given in Theorem 11 as 
Equation: 


CO 


So(- 5(Pcry)) = S> exp Gazu 
r=1 r=1 ie 
“il 
= dae 


T= 


= Ww, 


and therefore the non-homogeneous MC defined by { Por) a is weakly ergodic. Now we can 


use Theorem 12 to show that the MC is strongly ergodic by cee that 
Equation: 


S° ll *(r) — Mera) I], < c- 
r=1 


Since we know from earlier in the proof that 77) (w”) is increasing for w" € # and eventually 


decreasing for w" € #°, there exists a9 € N such that for any r; > ro, 
Equation: 


So il ty — meray I, = SS YS (ren (w”) - 7 (w”)) 
T=P0 we SE T=T0 


+ SY De (w") = ter41) (w")) 


w'¢H T=10 

= SO (mr41) (w") — mr) (w")) 
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+ S2 (mr) (w") — mer41) (w")) 
w' EH 
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Since the right hand side does not depend on 71, then we have that 
Equation: 


oe) 


| T(r) — W(r41) ll < 00. 
fol 


This implies that the non-homogeneous MC used by MCMC algorithm is strongly ergodic, and 
thus completes the proof of Theorem 10. 


Compression of nonparametric sources 


Consider a stationary ergodic source, where we no longer assume that it is parametric. We need to 
define notions of probability that fit the universal framework. For this, we study Kac's lemma [link]. 


We have a stationary source, {..., Xn, X-ni1,---,Xo0,X1,---}; 


z= aa = {X_941, X_s42,..., Xo}. Define N, as the number of shifts forward of a window of 
length £ until we see z again; this is called recurrence time. Define 
Equation: 
_ Jl Xp =2 
k= ’ 
0 else 


e.g., Yo = 1, then N, is the smallest positive number for which Y; = 1. Note that iy is a binary 
stationary ergodic source. Define Q, =Pr {Y, = 1;1 <j <k—1,Y; = 0|Yo = 1}. Then the 
average recurrence time can be computed, w = er £Q: = E|[N,]. We can now present Kac's 
lemma. 


Lemma 3 [link] = PAY=I) and E[N, = PP Gan — z| = Fay 

Let A = {Y,, = 1 for some — co < n < +oo}. Because z just appeared, then its probability is 
sic3 . __ Pr(A) 

positive. We will prove that = Pr(Yo=t)* 


Define Bt = {Y, =1,for some 0 <n < oo}, andB = {Y, = 1,for some n < 0}. Then 
A=BtUB =(BtNB )U(BtN(B-)°) U((B*)° NB). we claim that 
Pr (ee al (B-)°) =Pr ((B*)° ‘al B) = 0. This can be shown formally, but is easily seen by 


realizing that if z appears at any time n (say, positive) then it must appear at some negative time n 
with probability 1, because its probability is positive. 


Therefore, we have 
Equation: 


So Pr (¥,=1,Y.=1LY, =0,-k<n<j) 
k= 


> Pr (p= 1) Pra. Y, S04 << gy y= 1) 
k= 


j=0 k=1 
=r (Yo _ 1) S0 iQ: 
i=1 
=PriYyy= De 
Pr(A) 


Therefore, ps = . We conclude the proof by noting that Pr (A) = 0. 


Pr(Yo=1) 
Let us now develop a universal coding technique for stationary sources. Recall 

Ay = +E [- log (Pr (X a3) . The asymptotic equipartition property (AEP) of information 
theory [link] gives 

Equation: 


Pr (x! -F log (Pr (X{)) — He 


> s) < €(6, £), 


where limy_,.. € (6, 2) = 0. Define in this context a typical set T'(d, €) that satisfies 
Equation: 


g— (Het 6)e <Pr (Xf) = 9-H) 


For a typical sequence z, F [N,.| x! = z| < 2He+9), Then 
Equation: 


N,) 


Pr (“ere Ss 25) =Pr (2 €T (6,0) Pr Ge. 
log (N,) 
g 


> Hy +26|z € T (6, 0) 


+ Pr (z¢T(6,2)) Pr ( > Hy + 262 ¢ 7 (6,0)) 


git. AEP 
= e(d,£) — 0. 


Consider our situation, we have a source with memory of length n and want to transmit X ; 


1. Choose n = 2/(#e+28) 

2. For z= X. . find the value of N, if it appears in memory. 

3. If it appears, then transmit a 0 flag bit followed by the value of N,. 
4. Else transmit a 1 followed by the uncompressed z. 


Transmitting z via N, requires {log (V,.)] bits, and so the expected coding length is 
Equation: 


log (N;) 
Pr ( ? 
< (2+(H,+6))+2-" (1+ £ log a) + € (£,5) (2 + £ log a). 


< Hy 26) (1 + 1+ log (N,) + 2° (1 +£ log @)) + €(£,6)(1+1+ £log a) 


After we normalize by @, the per symbol length converges to 4 + Hy + 26. 


Note that this analysis assumes that the entropy Hy for a block of £ symbols is known. If not, then we 
can have several sets that are each adjusted for some different entropy level. 


Universal Coding of the Integers 


Coding of an index also appears in other information theoretic problems such as indexing a coset in 
distributed source coding and a codeword in lossy coding. What are good ways to encode the index? 


If the index n lies in the range n € {1,...,.N} and N is known, then we can encode the index using 
[log (V)]| bits. However, sometimes the index can be unbounded, or there could be an (unknown) 
distribution that biases us strongly toward smaller indices. Therefore, we are interested in a universal 
encoding of the integers such that each n is encoded using roughly log (7) bits. 


To do so, we outline an approach by Elias [link]. Let m be a natural number, and let b(m)be the binary 
form for n. For example, b(0) = 0, 6(1) = 1, (2) = 10, 6(3) = 11, and so on. Another challenge is 
that the length of this representation, |b(n)|, is unknown. Let u(k) = 0*11, for example 

u(4) = 0001. We can now define e(n) = u(|b(|b(m)|)|)b(/b(m) |)b(m). For example, for n=5, we have 
b(n) = 101. Therefore, |b(n)| = 3, b(|b(n)|) = 11, |b(/b(m)|)| = 2, u([b(|b(m)|)|) = 01, and 

e(n) = 0111101. We note in passing that the first 2 bits of e(n) describe the length of b(|b()|), 
which is 2, the middle 2 bits describe b(|b(m)|), which is 3, and the last 3 bits describe b(n), giving the 
actual number 5. It is easily seen that 

Equation: 


le(n)| = 2/6(/b(m)|)| + [6(m)| 
< 2 log [log (n+ 1) + 1] + log[n +1]. 


Sliding Window Universal Coding 


Having described an efficient way to encode an index, in particular a large one, we can employ this 
technique in sliding window coding as developed by Wyner and Ziv [link]. 


Take a fixed window of length n., and a source that generates characters. Define the history as x7", 
this is the portion of the data x that we have already processed and thus know. Our goal is to encode 


the phrase y; = ra ge where Ly is not fixed. The length L, is chosen such that there exists 


Equation: 
Y= On = Sam 
for some m € {0,1,...,n. — 1}. For example, if « = abcbababc and we have processed the first 5 


symbols abcba, then L; = 3 and m; = 1, where abcba is the history and we begin encoding after that 
part. 


The actual encoding of the phrase y; sends | f(y1)| bits, where 
Equation: 


_ f [log (nw)] + le (Z1)|, log (n2) < Li log (a) 
If ui )i= { [Ly log (a) | + |e (Z1)|, else 


First we encode Lj; then either the index m into history is encoded or y; is encoded uncompressed; 
finally, the first Z, characters ai ‘ are removed from the history database, and y; is appended to its 


end. The latter process gives this algorithm a sliding window flavor. 


This rather simple algorithm turns out to be asymptotically optimal in the limit of large n,,. That is, it 
achieves the entropy rate. To see this, let us compute the compression factor. Consider a” where 


N >> ny, R(N) = Be) | N= + 3 ly;|, and |y;|= L1, where C’ is the number of phrases 
required to encode ay . The number of bits / CaM ) needed to encoded ai using the sliding window 
algorithm is 

Equation: 


Cc 
L(y’) = [nw log (a)] + Sle (LZi)| + 1+ min {flog (n..)], [Li log (a)]}. 


i=1 


Therefore, the compression factor R (N) satisfies 
Equation: 


C 


= [Mw log (a) | ! = Done (nw)| +1 +11 log (Li), 72 log (Li)}- 


R(N) 7 


Claim 3 [link] As lim,,,-,. and limjy_;0, the compression factor R (NV) converges to the entropy 
rate H. 


A direct proof is complicated, because the location of phrases is a random variable, making a detailed 
analysis complicated. Let us try to simplify the problem. 


History Message N 


Block partitioning in analysis of sliding window 
algorithm [link]. 


— NywtN log(nw) 
Take lp that divides n., + N to create a H2e 


window size in the algorithm from [link]. This partition into blocks appears in [link]. Denote by 
4s 
Si —Io+1* 


intervals, where Ig = , which is related to the 


wi, (x) > 0 the smallest value for which x = x Using Kac's lemma [link], 
Equation: 


1 


n Pr (2: = z) 


Pr (wi, (a) > nei = z) < 


Claim 4 [link] We have the following convergence as / increases, 
Equation: 


l00 
Pr (wi (x) > a) a i, 


Because of the AEP [link], 
Equation: 


Therefore, for a typical input Pr (x!) > 2Q-lo(H+e) 


Recall that the interval length is lp = toa) , and so the probability that an interval cannot be found in 
the history is 
Equation: 
9—o(H+2€) i 
— 9—40€ 
9—lo(H+e) =2 70 


For a long enough interval, this probability goes to zero. 


Redundancy of parsing schemes 


There are many Ziv-Lempel style parsing algorithms [link], [link], [link], and each of the variants has 
different details, but the key idea is to find the longest match in a window of length n,,. The length of 


the match is L, where we remind the reader that LD ~ seer 


Now, encoding L requires log (n,,) bits, and so the per-symbol compression ratio is eel which in 
the limit of large n,, approaches the entropy rate H. 


However, the encoding of Z must also desribe its length, and often the symbol that follows the match. 
These require length log (ZL) ~log (log (n,)), and the normalized (per-symbol) cost is 
Equation: 


log (log (nw)) _ _,/( log (log (nw)) 
L 7 of log (nw) ) 


Therefore, the redundancy of Ziv-Lempel style compression algorithms is proportional to 


O (texte , which is much greater than the O (22s ) that we have seen for parametric sources. 


The fundamental reason why the redundancy is greater is that the class of non-parametric sources is 
much richer. Detailed redundancy analyses appear in a series of papers by Savari (c.f. [link]). 


Parsing for Lossy Compression 


The parsing schemes that we have seen can also be adapted to lossy compression. Let us describe 
several approaches along these lines. 

Fixed length: The first scheme, due to Gupta et al. [link], constructs a codebook of size ~ QL R(D) 
codewords, where L is the length of the phrase being matched and R(D) is the rate distortion 
function. The algorithm cannot search for perfect matches of the phrase, because this is lossy 
compression. Instead, it seeks the codeword that matches our input phrase most closely. It turns out 
that for large L the expected distortion of the lossy match will be approximately D per symbol. 


Variable length: Another approach, due to Gioran and Kontoyiannis [link], constructs a single long 
database string, and searches for the longest match whose distortion w.r.t. the input is approximately D 
; the location and length of the approximate match are encoded. Seeing that the database is of length 

~ 24 encoding the location requires ~ LR bits, and the D-match (a match with distortion D w.xr.t. 
the input string) is typically of length ~ L, giving a per-symbol rate of ~ R(D) bits. 


An advantage of the latter scheme by Gioran and Kontoyiannis [link] is reduced memory use. The 
database is a string of length ~ 2"R(D) instead of a codebook comprised of ~ 2"R(Y) codewords, 
each of length Z. On the other hand, the Gupta et al. algorithm [link] has better R.D performance, 
because it does not need to spend ~log (L) bits per phrase to describe its length. An improved 
algorithm, dubbed the hybrid algorithm by Gioran and Kontoyiannis, constructs a single database and 
performs fixed length coding for the best match of length L in the database. Therefore, it combines the 
memory usage of a single database approach with the RD performance of fixed length coding. 


Notation 


e x - input sequence 
e n- length of x 
‘2S Fito tee Seep d; 
e a - discrete alphabet 
e r- cardinality of a 
° mz (a) - the number of times that a € a appears in x 
e P, (a) - empirical probability 
e Q - the true i.i.d. distribution of x 
e H - entropy 
e D(- || -) - divergence 
e |-| - rounding up 
e R- rate of source code 
¢ Tg (6) - set of inputs that are 6-typical with respect to (w.r-t.) Q 
° @ -acode 
e (-)© - complement of set 
¢ -” _ transpose of vector or matrix 
e I(x) - length of code in bits 
e x’ = S,-stepVn € Z, 2), =2n11 
e @- parameters of parametric source (can contain multiple scalar 
parameters) 
e M - number of states of unifilar source 
e S = {1,...,M} - states of unifilar source 
© S1,..., $n, - State sequence 
e q(s|s’) - transition probability 
e St = g(St_-1, X¢-1) - next state function 
e A -class of parametric models 
* cy, - collection of lossless codes for length-n inputs 
e J(-;-) - mutual information 
rn (1, 0) - how far a coding length function J is from the entropy 
e h2(-) - binary entropy 
I(-) - Fisher information 
e J(-) - Jeffreys' prior 
e @yzr - Maximum likelihood parameter 
e R,, R* - min-max and max-min redundancy 


T - set of leaves of a context tree source 

$ - State of context tree 

Nz (s,a) - number of times that a € @ appears in context s in x 
Es optimal context tree source 

D - maximal depth of tree source 


ie - optimal tree structure to encode symbols whose context ends with 
8 

MDL(s) - minimum description length required for encoding these 
symbols 

KT(.,-) - Krichevsky-Trofimov coding length 

(y, 7) - BWT output consisting the permuted sequence and index 

x - approximate version of x 

q@ - reproduction alphabet 

d(-,-) - distortion metric 

d «x,x - distortion between sequences 

D - distortion level 

R(D) - minimal rate such that x can be described 

t - temperature in simulated annealing 

Z; - normalization factor in Boltzmann pmf 

N,. - recurrence time 


